/*
 * com.javable.cese -- Cell Electrophysiology Simulation Environment (CESE)
 *
 * Copyright (C) 1999 - 2005 Sergey Missan
 * 
 * Dalhousie Univ., Halifax, NS, Canada
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public
 * License along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 *
 */

package models.com.javable.cese.models;

import java.io.Serializable;

import com.javable.cese.templates.AbstractAPModel;

/**
 * The Hodgkin-Huxley Squid Axon Model, 1952. Detailed list of equations and
 * model description are provided in
 * <p>
 * J Physiol 1952;117:500-544.
 */

public class HodgkinHuxleyModel_1952 extends AbstractAPModel implements
        Serializable {

    /*
     * The beginning here lists all the variables and defines them (this code
     * uses all global variables)
     */

    /* Cell Geometry */
    private double l = 0.01; // Length of the cell (cm)

    public synchronized void setl(double sl) {
        l = sl;
        getCellGeometry();
    }

    public double getl() {
        return l;
    }

    private double a = 0.0011; // Radius of the cell (cm)

    public synchronized void seta(double sa) {
        a = sa;
        getCellGeometry();
    }

    public double geta() {
        return a;
    }

    private double vcell; // Cell volume (uL)

    public double getvcell() {
        return vcell;
    }

    private double ageo; // Geometric membrane area (cm^2)

    public double getageo() {
        return ageo;
    }

    private double acap; // Capacitive membrane area (cm^2)

    public double getacap() {
        return acap;
    }

    private double vmyo; // Myoplasm volume (uL)

    public double getvmyo() {
        return vmyo;
    }

    private double vmito; // Mitochondria volume (uL)

    public double getvmito() {
        return vmito;
    }

    private double vcleft; // Cleft volume (uL)

    public double getvcleft() {
        return vcleft;
    }

    /* Terms for Solution of Conductance and Reversal Potential */
    public static final double R = 8314; // Universal Gas Constant (J/kmol*K)
    public static final double frdy = 96485; // Faraday's Constant (C/mol)

    private double temp = 310; // Temperature (K)

    public synchronized void settemp(double stemp) {
        temp = stemp;
    }

    public double gettemp() {
        return temp;
    }

    /* Ion Valences */
    public static final double zna = 1; // Na valence
    public static final double zk = 1; // K valence

    /* Ion Concentrations */

    /* Myoplasmic Na Ion Concentration Changes */
    private double naiont; // Total Na Ion Flow (mM/ms)

    public double getnaiont() {
        return naiont;
    }

    private double dnai; // Change in Intracellular Na Concentration (mM)

    public double getdnai() {
        return dnai;
    }

    /* Myoplasmic K Ion Concentration Changes */
    private double kiont; // Total K Ion Flow (mM/ms)

    public double getkiont() {
        return kiont;
    }

    private double dki; // Change in Intracellular K Concentration (mM)

    public double getdki() {
        return dki;
    }

    private double er = -75.0; // Reversal Potential (mV)

    public synchronized void seter(double ser) {
        er = ser;
    }

    public double geter() {
        return er;
    }

    /* Fast Sodium Current (time dependent) */
    private double ina; // Fast Na Current (uA/uF)

    public synchronized void setina(double sina) {
        ina = sina;
    }

    public double getina() {
        return ina;
    }

    private double gna = 120.0; // Max. Conductance of the Na Channel (mS/uF)

    public synchronized void setgna(double sgna) {
        gna = sgna;
    }

    public double getgna() {
        return gna;
    }

    private double ena; // Reversal Potential of Na (mV)

    public double getena() {
        return ena;
    }

    private double ah; // Na alpha-h rate constant (ms^-1)

    public double getah() {
        return ah;
    }

    private double bh; // Na beta-h rate constant (ms^-1)

    public double getbh() {
        return bh;
    }

    private double am; // Na alpha-m rate constant (ms^-1)

    public double getam() {
        return am;
    }

    private double bm; // Na beta-m rate constant (ms^-1)

    public double gebm() {
        return bm;
    }

    private double h = 0.6; // Na inactivation h

    public double geth() {
        return h;
    }

    private double m = 0.05; // Na activation m

    public double getm() {
        return m;
    }

    private double scgna = 1;

    // Scaling factor for Max. Conductance of the Na Channel
    public synchronized void setscgna(double sscgna) {
        scgna = sscgna;
    }

    public double getscgna() {
        return scgna;
    }

    /* Time-dependent Potassium Current */
    private double ikt; // Time-dependent K Current (uA/uF)

    public synchronized void setikt(double sikt) {
        ikt = sikt;
    }

    public double getikt() {
        return ikt;
    }

    private double gkt = 36.0;

    // Channel Conductance of Time-dependent K Current (mS/uF)
    public synchronized void setgkt(double sgkt) {
        gkt = sgkt;
    }

    public double getgkt() {
        return gkt;
    }

    private double ekt; // Reversal Potential of K (mV)

    public double getekt() {
        return ekt;
    }

    private double an; // K alpha-n rate constant (ms^-1)

    public double getan() {
        return an;
    }

    private double bn; // K beta-n rate constant (ms^-1)

    public double getbn() {
        return bn;
    }

    private double n = 0.325; // Time-dependent K activation

    public double getn() {
        return n;
    }

    private double scgkt = 1;

    // Scale factor for Channel Conductance of Time-dependent K Current
    public synchronized void setscgkt(double sscgkt) {
        scgkt = sscgkt;
    }

    public double getscgkt() {
        return scgkt;
    }

    /* Leakage Current */
    private double il; // Leakage current (uA/uF)

    public synchronized void setil(double sil) {
        il = sil;
    }

    public double getil() {
        return il;
    }

    private double gl = 0.3;

    // Channel Conductance of Leakage Current (mS/uF)
    public synchronized void setgl(double sgl) {
        gl = sgl;
    }

    public double getgl() {
        return gl;
    }

    private double el; // Reversal Potential of Leakage (mV)

    public double getel() {
        return el;
    }

    private double scgl = 1;

    // Scaling factor for Channel Conductance of Leakage Current
    public synchronized void setscgl(double sscgl) {
        scgl = sscgl;
    }

    public double getscgl() {
        return scgl;
    }

    /* ** CONSTRUCTOR ** */
    public HodgkinHuxleyModel_1952() {
        super();
        getCellGeometry();
    }

    public void setCurrents() {
        currentINa();
        currentIKt();
        currentIL();
        currentIt();
    }

    public void doSimulate() {
        setCurrents();

        vnew = v - it * dt;
        dvdtnew = (vnew - v) / dt;
        t = t + dt;
        v = vnew;
        dvdt = dvdtnew;
    }

    /** ***************************************************** */
    /**
     * Cell Geometry
     */
    protected void getCellGeometry() {
        vcell = 1000 * Math.PI * a * a * l;
        ageo = 2 * Math.PI * a * a + 2 * Math.PI * a * l;
        acap = ageo * 2;
        vmyo = vcell * 0.68;
        vmito = vcell * 0.26;
        vcleft = vcell * 0.12 / 0.88;

        // set initial variable values
        setetime(100.0);
        setdt(0.01);
        v = -68.0;
        bcl = 50.0;
    }

    /* Functions that describe the currents begin here */

    /**
     * Calculates Fast Na Current
     */
    protected void currentINa() {
        gna = 120.0 * scgna;

        if (v == -35.0) {
            am = 1.0;
        } else {
            am = (0.1 * (v + 35.0)) / (1.0 - Math.exp((-0.1 * (v + 35.0))));
        }

        bm = 4.0 * Math.exp((v + 60.0) / -18.0);

        ah = 0.07 * Math.exp((v + 60.0) / -20.0);
        bh = 1.0 / (Math.exp((-0.1 * (v + 30.0))) + 1.0);

        if (t == 0) // Steady State Values
        {
            h = ah / (ah + bh);
            m = am / (am + bm);
        } else {
            h = ah / (ah + bh) - ((ah / (ah + bh)) - h)
                    * Math.exp(-dt / (1 / (ah + bh)));
            m = am / (am + bm) - ((am / (am + bm)) - m)
                    * Math.exp(-dt / (1 / (am + bm)));
        }

        ena = er + 115.0;
        ina = (gna * m * m * m * h) * (v - ena);
    }

    /**
     * Calculates Time-dependent K Current
     */
    protected void currentIKt() {
        gkt = 36.0 * scgkt;

        if (v == -50.0) {
            an = 0.1;
        } else {
            an = (0.01 * (v + 50.0)) / (1.0 - Math.exp((-0.1 * (v + 50.0))));
        }

        bn = 0.125 * Math.exp((v + 60.0) / -80.0);

        if (t == 0) // Steady State Values
        {
            n = an / (an + bn);
        } else {
            n = an / (an + bn) - ((an / (an + bn)) - n)
                    * Math.exp(-dt / (1 / (an + bn)));
        }

        ekt = er - 12.0;
        ikt = gkt * Math.pow(n, 4.0) * (v - ekt);
    }

    /**
     * Calculates Leakage Current
     */
    protected void currentIL() {
        gl = 0.3 * scgl;
        el = er + 10.613;
        il = gl * (v - el);
    }

    /**
     * Calculates Total Current
     * 
     * Total sum of currents is calculated here, if the time is between stimtime =
     * 0 and stimtime = tstim, a stimulus is applied
     */
    public void currentIt() {
        naiont = ina;
        kiont = ikt + il;

        if ((t >= tstim) && (t < (tstim + dt))) {
            stimtime = 0;
            tstim = tstim + bcl;
            dvdtmax = false;
        }

        if ((stimtime >= 0) && (stimtime < stl)) {
            it = st + naiont + kiont;
        } else {
            it = naiont + kiont;
        }

        stimtime = stimtime + dt;
    }

}
